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Abstract 

In this paper, we investigate the motion of a neutrally buoyant cylinder of circular 
or elliptic shape in two dimensional shear flow of a Newtonian fluid by direct nu- 
merical simulation. The numerical results are validated by comparisons with existing 
theoretical, experimental and numerical results, including a power law of the normal- 
ized angular speed versus the particle Reynolds number. The centerline between two 
walls is an expected equilibrium position of the cylinder mass center in shear flow. 
When placing the particle away from the centerline initially, it migrates toward an- 
other equilibrium position for higher Reynolds numbers due to the interplay between 
the slip velocity, the Magnus force, and the wall repulsion force. 

keywords: Shear flow; Neutrally buoyant particle; Equilibrium height; Fictitious do- 
main/distributed Lagrange multiplier method; Finite element. 



1. Introduction 

The problem of particle motions in shear flows is crucially important in many 
engineering fields such as the handling of a fluid-solid mixture in slurry, colloid, and 
fluidized bed. Segre and Silberberg J27J |2B] experimentally studied the migration of 
dilute suspensions of neutrally buoyant spheres in a tube Poiseuille flow. The particles 
migrate away from the wall and the centerline to accumulate at about 0.6 of the tube 
radius from the centerline. The experiments of Segre and Silberberg [27J EE] have 
had a large influence on fluid mechanics studies of migration and lift of particles. 
Comprehensive reviews of experimental and theoretical works have been given by 
Brenner [2], Cox and Mason [5J, Feuillebois [TT] and Leal [18] . 
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Among the theoretical studies of the neutrally buoyant particle migration in linear 
shear flow, Bretherton [3J found an expression for the lift force per unit length on a 
cylinder in an unbounded two-dimensional linear shear flow at small Reynolds num- 
ber. Saffman's lift force [26] on a sphere of radius a in an unbounded linear shear 
flow with shear rate G is F s = 6A6pVa 2 (Gu) 1 ^ 2 = 6A6puaV(Re P Y^ 2 where v is the 
kinetic viscosity of the fluid, p is the density of the fluid, Re p = Ga 2 / v is the particle 
Renolds number, and V is the slip velocity of the sphere. In a bounded linear shear 
flow, Ho and Leal [16] examined the motion of a rigid sphere with inclusion of the 
inertia effects at small Reynolds numbers by a regular perturbation method. The 
sphere reaches a stable lateral equilibrium position which is the midway between the 
walls. Vasseur and Cox [22] also obtained the same stable lateral equilibrium posi- 
tion. Ho and Leal require that Re p /n 2 <C 1 which is more restrictive than the one 
Re p / k <C 1 required by Vasseur and Cox where k = 2a /H is the confined ratio, H 
being the distance between two walls. Direct numerical simulations have been used 
for understanding particle motion in shear flows. Feng et al. [10] investigated the mo- 
tion of neutrally buoyant and non-neutrally buoyant circular particle in plane shear 
and Poiseuille flows using a finite element method and obtained qualitative agreement 
with the results of perturbation theories and of experiments. The numerical results 
of a neutrally buoyant circular cylinder in a shear flow of Re p = 0.625 have been dis- 
cussed in details. The cylinder migrates back to the midway between two walls due 
to the wall repulsion at the small Reynolds number. They have suggested that that 
three factors, namely the wall repulsion due to a lubrication effect, the slip velocity, 
and the Magnus type of lift, are possible responsible for the lateral migration. Ding 
and Aidun [8] studied numerically the dynamics of a cylinder of circular or elliptic 
shape suspended in shear flow at various particle Reynolds number. They obtained 
the transient from being rotary to stationary as the particle Reynolds number is in- 
creased for an elliptic cylinder. For the cases of the circular cylinder, the effect of the 
two walls on the rotation speed u has been studied. For the confined ratio k = 0.5, 
Ding and Aidun reported \w\/G oc Re p - 28 . Similar result, \co\/G oc Re p 0,25 , has been 
observed experimentally by Zettner and Yoda [31] . 

In this paper, we first discuss the generalization of a distributed Lagrange multi- 
plier/fictitious domain method (DLM/FD method) developed in [22J to non-spherical 
neutrally buoyant cylinders in two-dimensional shear flows and to the cases with zero 
angular velocity as a constraint. Similar DLM/FD methods for freely moving neu- 
trally buoyant particle has been successfully applied, in [201 E3] |3Q] , to simulate par- 
ticulate flow in two and three dimensions with neutrally buoyant particles. But for 
the cases of a neutrally buoyant elliptic cylinder in two dimensional flows, the method- 
ology has not been fully validated yet. We have validated the numerical method by 
comparing with the computational results in Ding and Aidun [S] for a cylinder of 
circular and elliptic shape and the experimental results in Zettner and Yoda [31] for 
a circular cylinder. On the wall effect on the angular velocity of the circular cylin- 
der, we have obtained oc Re~ ' 2771 for the confined ratio k = 0.5 which is in 
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a good agreement with the results obtained by Ding and Aidun [5] and Zettner and 
Yoda [31] . We have also studied the wall effect on the angular velocity of the elliptic 
cylinder which is more complicated due to the non-circular shape. 

Concerning the equilibrium position of a neutrally buoyant circular cylinder in shear 
flow, recent studies by Cherukat, McLaughlin and Dandy j5] and Kurose and Komori 
[TT] focus on lift and drag on a stationary sphere in unbounded linear shear flow. The 
equilibrium positions have not been studied in these works. Feng and Michaelides [H] 
have investigated the equilibrium positions of non-neutrally buoyant circular cylinders 
in two-dimensional shear flow. In their simulations, the density ratio between the solid 
and fluid is between 1.005 and 1.1. The equilibrium heights of their lightest circular 
cylinder (the density ratio of 1.005) are far below the centerline. We have obtained 
that the cylinder stays in the middle between two walls as expected when placing it 
there initially. But when the initial position of the mass center of a circular cylinder is 
away from the centerline, the equilibrium position depends on the particle Reynolds 
number Re p and the confined ratio k. For those placing away from the centerline 
initially, the circular cylinder migrates back to the centerline for Re p < Re P)C where 
Re PtC is the critical value. For Re p > Re p>c , the equilibrium position is between the 
wall and the centerline. The critical particle Reynolds number is increasing when 
increasing the confined ratio. Concerning the Magnus lift effect on the equilibrium 
position, we have added a constraint, zero angular velocity, to the motion of a circular 
cylinder and obtained that the equilibrium position of the circular cylinder moving 
with zero angular velocity is lower than those of freely moving cylinder when both are 
away from the middle. These results show that the Magnus lift does play a role as 
expected. Also from the computed slip velocity of the circular cylinder, it shows that 
the circular cylinder lags the fluid, at least for Re p > Re PjC , which means that there is 
a force pushing the cylinder toward the wall (see Fig. [Ijfor the setup of the boundary 
conditions). Hence the equilibrium position of the cylinder is up to the interplay 
between the slip velocity, the Magnus lift and the wall repulsion force. The content 
of this paper is as follows: In Section 2 we introduce a fictitious domain formulations 
of the model problem associated with the neutrally buoyant long particle cases; then 
in Section 3 we discuss the time and space discretization and in Section 4 we present 
and discuss the numerical results. 



2. A fictitious domain formulation of the model problem 

A fictitious domain formulation with distributed Lagrange multipliers for flow 
around freely moving particles and its associated computational methods have been 
developed and tested in, e.g., [T21 [331 EU 125]- For the cases of neutrally buoyant 
particles, similar methodologies have been developed in [20l E21 |23l [30]. But for the 
cases of a neutrally buoyant elliptic cylinder in two dimensional flows, the methodol- 
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Figure 1: An example of two-dimensional flow region with a rigid body. 

ogy has not been discussed and fully validated yet. In this paper, we first discuss the 
formulation for the case of a neutrally buoyant elliptic cylinder and then present the 
numerical results to validate the methodology. 

Let f2 C 1R 2 be a rectangular region filled with a Newtonian viscous incompressible 
fluid (of density p and dynamic viscosity p) and containing a freely moving neutrally 
buoyant rigid particle B centered at G = {G±, G^}* of density p. The flow is modeled 
by the Navier-Stokes equations and the motion of the particle B is described by the 
Euler-Newton's equations. We define 

W S0:P = {v|v G (if 1 (f2)) 2 , v = g on the top and bottom of and 

v is periodic in the x\ direction}, 
W 0jP = {v|v G (if^fi)) 2 , v = on the top and bottom of and 

v is periodic in the X\ direction), 

L 2 = {q\qeL 2 (Q), [ grfx = 0,}, 
Jn 

A (t) = {ii\ii G (H\B(t))) 2 , < it, e, > B(t) = 0, % = 1, 2, < fi, > B(t) = 0} 

with ei = {1, 0}*, e 2 = {0, 1}*, Gx = {— (x 2 — G 2 ), X\ — G{\ 1 and < ■, • >B(t) & n inner 
product on Ao(t) which can be the standard inner product on (if 1 (i?(t))) 2 . For simple 
shear flow, we have go = (U /2, 0)* on the top wall and (—U/2, 0)* on the bottom wall. 
Then the fictitious domain formulation with distributed Lagrange multipliers for flow 
around a freely moving neutrally buoyant particle of the elliptic shape is as follows 



For a.e. t > 0, find u(t) G W So , p , p{t) G L%, V G (t) G IR 2 , G(f) G H 2 , 
w(t) G IR, 6{t) G IR, \{t) G A (t) such that 



T.-W. Pan & et al. 



(2) 
(3) 
(4) 

(5) 
(6) 

(7) 



du 



+ (u- V)u 



v <ix + fi / Vu : Vv dx. — I pV • v rfx 

<n Jn 



=< A,v > B(t) , Vv G W ,p, 
qV ■ u(t)dx = 0, Vg G L 2 (fi), 



V ( 



< j*,u(i) > B (t)= 0, Vjz g Aq(*), 
dG 

V G (0) = V° G , w(0) = a; , G(0) = G° = {G°, G°}<, 0(0) = 0°, 



u(x,0) = u (x) 



u (x), Vx G Q\B(0), 

V G + w°{-(x 2 - G° 2 ), X! - G?}*, Vx G 5(0), 



where u and p denote velocity and pressure, respectively, A is a Lagrange multiplier, 
Vg is the translation velocity of the particle B, u is the angular velocity of B, and 8 
is the angle between the horizontal direction and the long axis of the elliptic cylinder 
sec Fig|TJ). We suppose that the no-slip condition holds on dB. We also use, if 
necessary, the notation (f>(t) for the function x — > <p(pc,t). 

Remark 1 . The hydrodynamical forces and torque imposed on the rigid body by the 
fluid are built in @-((7]) implicitly (see [HI H3] for details), thus we do not need to 
compute them explicitly in the simulation. Since in ([l])-([7]) the flow field is defined 
on the entire domain Q, it can be computed with a simple structured grid. 

Remark 2. In (|3]), the rigid body motion in the region occupied by the particle is 
enforced via the Lagrange multiplier A. To recover the translation velocity V G (t) 
and the angular velocity u(t), we solve the following equations 



(9) 



< e h u(t) - V G (t) - u(t) Gx > B{t) = 0, for i = 1, 2, 

< Gx^, u(t) - V G (t) - w(t) Gx^ > B( t)= 0. 



Remark 3. To investigate the effect of the Magnus type of lift on the lateral migration 
of the cylinder, we have considered the cases of the cylinder freely moving in shear 
flow with zero angular velocity. For this special consideration, we have the following 
modified formulation 

For a.e. t > 0, find u(t) G W g0jP , p{t) G L 2 , V G (t) e IR 2 , G(t) G 1R 2 , 
\(t) G Ao(t) such that 
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(10) 

(11) 
(12) 

(13) 
(14) 

(15) 



du 



=< A,v > B{t) , Vv e W ,p, 
qV ■ u(t)dx = 0, Vg G L 2 (Q), 



v dx. + jx I Vu : Vv dx — / • v dx 



< M,u(fj > B (t)= 0, V/x G A (t), 

^ = V G 

dt G ' 

V G (0)=V G , G(0) = G° = {G?,G^, 



u(x, 0) = u (x) 



u (x), Vx G n\B(o), 



V° G , Vx G P(0) 
with the modified multiplier space 

A (t) = G (P 1 ^))) 2 , < M , e, > B{t) = 0,t = l, 2}. 

The translation velocity Vg(£) is recovered via 
(16) < e t , u(t) - V G (t) > B{t) = 0, for i = l, 2. 



3. Space approximation and time discretization 

Concerning the space approximation of problem (JI])- ^ by a finite element method, 
we have used P1-1S0-P2 and P\ finite elements for the velocity field and pressure, 
respectively (like in Bristeau et al. [lj). We approximate then W SOtP , Wq )P , L 2 and L\ 
by the following finite dimensional spaces 

(17W S0 ,h = K I v h G (C (H)) 2 , v h | T G Pi x P 1? VT G 7;, v h = g on the top 
and bottom of Q and v is periodic at T in the x± direction }, 

(18) WV = {v h I v,, G (C (H)) 2 , v fc | T e Pi x Pi, VT G T ft , v ft = on the top 
and bottom of Q and v is periodic at T in the x\ direction }, 



(19) 

and 
(20) 



{Qh 



q h 6 C°(n), q h \ T G Pi, VT G 7k, <?h «s periodic 
at T m £/ie £1 direction}, 



T 2 



{q h \qh G L£, / q h dx = 0}, 



respectively; in (17)-(20), Pi is the space of polynomials in two variables of degree 
< 1. 



T.-W. Pan & et al. 



7 



A finite dimensional space approximating A (t) is defined as follows: let {xj}^ be 
a set of points covering B(t); the discrete multiplier space A^(t) is defined by 



(21) 



A 



EN 



where £(■) is the Dirac measure at x = 0. Then, we have the inner product defined 
by 



(22) < fJ, h , v h >B h {t)= 5^=i^ ' Vft ( x ')' V ^ e V/l G 
and approximate A (t) by 

(23) A 0A (t) = {fi h \fi h e A h (t), < n h , ei > Bfl (t)= 0, £ = 1,2, </^,Gx >B h (t)= 0}. 

Using the above finite dimensional spaces leads to the following approximation of 
problem Q-0: 

For a.e. t > 0, /ind u(t) G W g0!ft , p(t) G L 2 0h , V G (t) G 1R 2 , G(t) G IR 2 , 
G IR, 9{t) G IR, A fc (*) G Ao,fc(t) swc/i &trf 



(24) 

(25) 
(26) 
(27) 

(28) 

(29) 
(30) 



at 



(u ft • V)u A 



v dx + ji / Vu fi : Vv c?x 



- / Ph^ ■ v dx =< A fe , v > Bh (t), Vv G W ,/i, 
/ qV ■ u h {t)dx = 0, Vg G L 2 



- 2 



< ji,u ft (t) > Bh (t)= 0, V/x G A , h (t), 

^=V G 

dt 

de 

~dt 



co, 



V Q (0) = V G , W (0) = co°, G(0) = G° = {G?, G°}*, (9(0) = 0°, 
u fc (x, 0) = Uo,fc(x) (withV • u 0A = 0). 



Applying a first order operator splitting scheme a la Marchuk-Yanenko [19] to the 



equations (24)- (30) at each time step and the Euler backward method in time for 



some subproblems, we obtain (after dropping some of the subscripts h): 



(31) 



u° = Uo,ft, V G , co°, G°, and 9° are given; 
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For n > 0, knowing u n , Vq, u n , G n , and 6 n , compute u n+1 / 6 and p n+l / 6 via the 
solution of 



(32) 



u n+l/6 _ u ra r 

PJ ^ vdx-|p" +1 / 6 V-vdx = 0, Vve%, 

/ gV • u" +1 / 6 rfx = 0, Vg G L*; u" +1 / 6 G W^,,, K +1 / 6 G L\ h . 



Then compute 

u n+2/6 ma 

the solution of 



(33) 
(34) 



fy^-vrfx + y (u" +1 / 6 - V)u-vrfx = 0, VvG%, on (t n ,t n+1 ), 

\ u(t n ) = u n+1 / (i - n u(t)eW S0:h: 



u 



n+2/6 _ „/>n+l 



= u(t n+1 ). 



Next, compute u ra+3//6 wa the solution of 
(35) 



u n+3/6 _ n+2/6 /• 

p / v dx + an / Vu n+3/6 • Vv dx = 0, 



Now predict the position and the orientation of the particle via: 

dG 



(36) 

(37) 
(38) 



dt 
d6 

~di 



V&/2, on (f\t n+1 ) 



w n /2, on (t n ,r +1 ) 



G(t n ) = G n , 0(t n ) = 9 n . 



Then set G n+4 / 6 = G(t n+1 ) and # n+4 / 6 = #(r +1 ). 

Now, compute u n+5 / 6 ; A n+5/6 7 V™ +5/6 7 and u n+5 / 6 via the solution of 



(39) { 



( f -.jrt+5/6 _ n n+3/6 r 

p I — vdx + /3u / Vu n+5 / 6 - V 

Jn 



n At 

< A, v > R n+4/6 , Vv G W oh , 



vdx 



^ ,, „n+5/6 ^ , — O Vf. c A™ +4 / 6 - ,,™+5/6 c w \n+5/6 c a "+4/6 
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and solve for Vq +5 ' 6 and c<; n+5//6 from 



(40) 



< ei,u n+5 / 6 - V 



n+5/6 



->-L 



5/6 G n+4 / 6 x > B n+4 /G = 0, for i = 1, 2, 



< G n+4//6 x u n+5//6 — V n+5,/6 — o; n+5 / 6 G n+4/,6 x > 



B. 



Finally, correct the position and the orientation of the particle via: 



(41) 

(42) 
(43) 

and set G 
V 



dG 

~dT 
d9 

dt ' 



V 



n+5/6 



/2, on (t n ,t n+1 ) 



w n+5/6 /2, on (t n ,t n+1 ) 



G(t n ) = G n+4/6 , 6(t n ) = 6 



n+4/6 



n+1 



n+5/6 



G(t 



n+l N 



an 



d 9 n+1 = 6(t n+1 ). Finally, we set u n+1 = u n+5 / 6 ; V G +1 



, andu n+l = w n + 5 / 6 . 



n+s 



In above algorithm Q-Q, we have t n+s = (n + s)At, A^+ s = A 0ih (t n+S ), B 
is the region occupied Dy the particle centered at G n+S . Finally, a and (3 verify 
a + /3 — 1; we have chosen a = 1 and /3 = in the numerical simulations discussed 
later. 



At each time step we have a sequence of subproblems (32), (33), (35) and (39). 



The degenerated quasi-Stokes problem (32) is solved by a preconditioned conjugate 
gradient method introduced in 



in which discrete elliptic problems from the pre- 
conditioning are solved by a matrix-free fast solver from FISHPAK by Adams et al. 



in pp. The advection problem (33) for the velocity field is solved by a wave-like equa- 



tion method as in [21]. Problem (35) is a classical discrete elliptic problem which 
can be solved by the same matrix-free fast solver. To enforce the rigid body motion 
inside the region occupied by the particles, we have applied the conjugate gradient 
method discussed in, e.g. 



Remark 4. To solve the problem ( 10 )-( 16 ) for the cases with zero angular velocity, we 



have an analogue algorithm by dropping 9 and u> from the algorithm (31)-(43) and 



replacing the equation (40) by 



< ei, u n+5/6 - V G +5/6 > B n+4 /6 = 0, for i = 1,2 



with the discrete multiplier space defined by 



(44) 



h 0> h(t) = {VhlVh e A h (t), < fi h ,ei >B h (t)= 0, % = 1,2}. 



4. Results and discussions 
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t 



Figure 2: Validation of the computational results by comparing both angle of incli- 
nation and angular velocity of an elliptical cylinder, (i) The angle 8/tt: (a) our result 
(solid line) at Re p = 0.02, (b) Ding's result (dash line) at Re p = 0.02, (c) our result 
at Re p = 0.25, (d) Ding's result at Re p = 0.25. (ii) The angular velocity u/ti: (e) our 
result at Re p = 0.02, (f) Ding's result at Re p = 0.02, (g) our result at Re p = 0.25, 
(h) Ding's result at Re p = 0.25. Jeffery's solutions at Re p = 0, Jq and J u , are also 
plotted for comparisons. 



4.1. The motion of an elliptic cylinder in linear shear flow 

We have first considered the motion of a neutrally buoyant elliptical cylinder in 
linear shear flow studied by Ding and Aidun [8] for the validation purpose. Their 
results were obtained by the lattice Boltzmann equation. The domain of computation 
is Q = [0, 5] x [0, 1], then the height is H — 1. The confined ratio is k = 2a/ H = 0.2, 
and the aspect ratio is AR = b/a = 0.5 where a is the length of the semi-major axis 
and b is the length of the semi-minor axis. The density of the fluid is p = 1 and 
the kinetic viscosity is determined by the specified value of the Reynolds number via 
v = Ga 2 /Re p . The shear rate is fixed at G = 1 and U = —GH = —1 (so the moving 
directions of the two walls are the same as those in Ding and Aidun |8J). The mesh 
size is h = 1/320 and the time step is At = 0.001. 

In Fig. [2j the computational results of the angle and angular velocity of the elliptic 
cylinder are in a good agreement with the Jeffery's solution [JS] for Re p = and the 
ones obtained by Ding and Aidun for Re p = 0.02 and 0.25, respectively. Fig. [3] shows 
that at various particle Reynolds numbers from to 7.5, our results of angular velocity 
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Figure 3: Angular velocity u/tt of an elliptical cylinder in shear flow at different Re p . 
(a) Our result (solid line) at Re p = 3.75; (b) Ding's result (dash line) at Re p = 3.75; 
(c) Our result at Re p = 7; (d) Ding's result at Re p = 7; (e) Our result at Re p = 7.5 
and (f) Ding's result at Re p = 7.5. Jeffery's solution J w at Re p = is also plotted for 
comparisons. 



match very well with the ones obtained by Ding and Aidun |Sj. The motion of the 
ellipse for Re p < 7 is a periodic rotation with non-uniform angular velocity, while for 
Re p > 7.5, the ellipse does not rotate at all; instead it takes a stationary orientation 
in shear flow (see, e.g., Ding and Aidun [8] for further details). An example of the 
velocity field of Re p = 8.25 around an elliptic cylinder with a stable orientation is 
shown in Fig. |4} 

The minimum angular velocity decreases as the particle Reynolds number is in- 
creased with a nearly straight line relationship as shown in Fig. [5j where Re p c = 
7.25 is the critical particle Reynolds number above which the rotational motion is 
stopped. Fig. [5] shows that the period of rotation increases rapidly as the particle 
Reynolds number is increased close to Re p ^ c , and the results are in good agreements 
with those obtained by Ding and Aidun [8]. The further studies of the motion of a 
cylinder of axisymmetric shape in simple shear flow will be reported in a forthcoming 
paper. 

4.2. The rotation of a cylinder in a simple shear flow 

In, e.g., Aidun and Ding [8J, the angular velocity of a circular cylinder suspended 
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Figure 4: The velocity field of the case of k—0.2 and Re p = 8.25. 




Figure 5: Minimum angular velocity u/n versus Re p of an elliptical cylinder (left). 
Period of motion T of an elliptical cylinder. It is noted that T increases to infinity as 
Re p approaches the critical value Re Pfi ~ 7.25 (right). 
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Figure 6: The velocity field of the case k=0.5 and Re p = 15. 
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Figure 7: Log-log plot of normalized angular velocity u/G vs Re p : (a). k=0.5 (blue 
*), 0.25 (blue x), 0.125 (blue □). (b). Ding & Aidun's results (2000): « = 0.5 (red 
+) and 0.25 (red x). (c). Zettner & Yoda's results (2001): k = 0.5 (red 0). 



in a simple shear flow has been studied via the direct numerical simulation. Experi- 
mental results have been also obtained in, e.g., Zettner and Yoda [31]. In this section, 
we focus on the wall effect of the rotation speed of a cylinder freely suspended in a 
simple shear flow with various confined ratios. We have first considered the cases of 
a circular cylinder suspended in the middle between two walls initially. The domain 
of computation is Q = [0, L] x [0, 1] with L = 16a where a is the circular cylinder 
radius. The density of the fluid is p = 1 and the kinetic viscosity is v = 0.012. The 
confined ratios are k = 2a/H = 2a =0.125, 0.25 and 0.5. We have varied the values 
of the shear rate G to have different values of the particle Reynolds number Re p . The 
initial position of the cylinder is at the midway between two walls. The mass center 
of the circular cylinder stays at the centerline between two walls in the simulations 
without giving any extra conditions for keeping it there. A typical velocity field (here, 
Re p = 15 and k=0.5) is shown in Fig. [6j At the zero Reynolds number, the rotation 
speed of a circular cylinder is G/2 from the Jeffery's solution. At the small particle 
Reynolds numbers, the ratio \oj\/G is about to converge to 0.5 when decreasing the 
confined ratio k from 0.5 to 0.125 as shown in Fig. [7j When k =0.5, the normalized 
angular speed \w\/G is found to be about 0.420 for Re p < 2. When k =0.25, \oj\/G is 
about 0.482 for Re p < 0.2. Both are very close to those obtained in Ding and Aidun 
[Sj. When k =0.125, \u\/G = 0.496 is much closer to 0.5 for Re p < 0.05 due to weaker 
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Figure 8: Log-log plot of normalized angular velocity \u m i n \/G versus Re p for a fixed 
confined ratio k = 0.2 (left) and a case of a fixed aspect ratio AR = 0.5 (right). 



effect from the walls. 

When increasing the particle Reynolds numbers, the log-log plot of \oj\/G versus 
Re p in Fig. M shows that \u\/G oc Re' - 2771 for 10 < Re p < 29 when « =0.5. The 
one obtainecfin Ding and Aidun jH] is \w\/G oc Re~ ' 28 and the experimental results 
reported in Zettner and Yoda [31] is \oj\/G oc Re~ ' 25 for the same confined ratio 
k =0.5. For the smaller values of the confined ratios, k =0.25 and 0.125, we have 
obtained \u)\/G oc i?e~ a2976 and \uj\/G oc -Re" 4310 respectively. The results of k =0.5 
and 0.25 are in a good agreement with those obtained in Ding and Aidun [S]. 

For an elliptic cylinder suspended in linear shear flow, we have studied the effect of 
the aspect ratio AR. The mass center of an elliptic cylinder also stays at the middle 
between two walls if it is placed there initially. Since the elliptic cylinder has zero 
angular velocity when the particle Reynolds number is larger than the critical value 
Re P}C , the behavior of the rotation of an elliptic cylinder is different from that of the 
circular cylinder. In Fig. [8j the minimal angular velocity decreases rapidly to about 
zero when the particle Reynolds number Re p is closer but less than Re PtC since the 
motion of the elliptic cylinder is about to transit into the one with a fixed orientation 
in linear shear flow. The minimal angular velocity decreases faster for the smaller 
aspect ratio shows that the slender shape is easier to reach a stable orientation in 
linear shear flow. But its maximal angular velocity shown in Fig. [8] have different 
behavior since the maximal angular velocity happens when the direction of the long 
axis is about perpendicular to the shear direction. For the small particle Reynolds 
numbers, the maximal and minimal values of the angular velocity are close to the 
Jeffery's solution as in Fig. 18] 
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Figure 9: The equilibrium height of the mass center of a circular cylinder versus Re p 
for k = 0.25. 

4.3. The equilibrium position 

In Ho and Leal [16] and Vasseur and Cox [29], they concluded that the sphere 
reaches a stable lateral equilibrium position which is the midway between the walls 
for small particle Reynolds numbers. In Feng et al. [10] . the cylinder migrates back 
to the midway between two walls at Re p = 0.625 when placing it away from the 
middle between two walls. Feng et al. have suggested that that three factors, namely 
the wall repulsion due to a lubrication effect, the slip velocity, and the Magnus type 
of lift, are possible responsible for the lateral migration. In the previous section, we 
have obtained that the centerline is always (at least in the range we have studied in 
this paper) the equilibrium position for the cylinder of elliptic or circular shape when 
it is positioned there initially. When placing the mass center initially away from the 
middle between two walls, the mass center may not migrate back to the centerline. 

For the cases of a circular cylinder of the confined ratio n = 0.25 as in the previous 
section, the radius of the circular cylinder is a = 0.125, the density of the fluid is 
p = 1, and the kinetic viscosity is v = 0.012. To check the effect of the length L of 
the channel, the initial height yinit-, the mesh size h and time step dt on the equilibrium 
height of the mass center, we have considered different sets of parameter values as 
indicated in Fig. [9} The final equilibrium heights in Fig. [9] are almost the same for 
each value of Re p considered here except those with zero angular velocity constraint. 
When Re p =1 and 2, the mass center of the freely moving cylinder migrates back to 
the middle between two walls. But for higher particle Reynolds numbers, we have 
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Figure 10: The equilibrium height of the mass center of a circular cylinder versus Re p 
for various values of k. 
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Figure 11: The comparison of the equilibrium height of the mass center of a non- 
neutrally buoyant circular cylinder and that of a neutrally buoyant circular cylinder 
versus Re p . 
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Figure 12: The histories of the slip velocity (top), the angular velocity (middle), and 
the height of the mass center (bottom) of the cylinder for Re p =1 and 5: free motion 
(solid blue line) and zero angular velocity (red dashed line). 



obtained an equilibrium height which is between the centerline and the wall. To find 
out the effect of the walls on the final equilibrium height, we have varied the distance 
H between two walls. The initial height yi n n is always 0.1 unit below the centerline. 
The length of the computational domain is L = 2. The equilibrium height versus 
the particle Reynolds number is shown in Fig. 10 As expected, the wall repulsion 



force is stronger for the larger confined ratio. Hence the cylinder migrates back to 
the midway between two walls for k = 0.5 for the range of Re p considered here and 
the critical particle Reynolds number Re PjC is increasing when increasing the confined 
ratio. We believe that the symmetric breaking shown in Figs. [9] and [10] is not a 
numerical artifact. To further validate this phenomenon, we have compared ours 
with the results obtained in Feng and Michaelides [9] where only the non-neutrally 




Figure 13: A snapshot of the velocity field of Re p = 5 obtained by following the 
cylinder mass center which is moving to the left: freely motion (top) and zero angular 
velocity (bottom). 



circular cylinders were considered. For the results presented in Fig. 11, the radius 
of circular cylinder is a = 0.05, the initial height is yi n n = 0.1, the kinetic viscosity 
is v — 0.05, and the length and the height of the computation domain are L = 2 
and H — 1, respectively. The density p s of the non-neutrally circular cylinder is 
either 1.002 or 1.005. We have applied the DLM/FD method developed in, e.g., 
[I2| [13] to obtain the computational results for the cases of the non-neutrally circular 
cylinder. Our results are in a good agreement with the results obtained in Feng and 
Michaelides [D] . The transition of the equilibrium height from the one associated with 
a non-neutrally buoyant cylinder to that of a neutrally buoyant cylinder is consistent 
and support the existence of the symmetric breaking of the equilibrium height. 

For the effect of the Magnus lift associated with the rotation, we have considered 
the cases of a circular cylinder of k =0.25 moving with zero angular velocity. The 
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equilibrium height of non rotating cylinder is lower than the one of the freely moving 
cylinder for Re p > 2 as in Fig. [9j These results indicate that the Magnus lift does play 
a significant role as expected. The similar results concerning the Magnus lift have 
also been observed in Feng and Michaelides [9]. When Re p —1, the wall repulsion 
force is strong enough to push the cylinder back to the centerline even without the 
help from the Magnus lift. 

We now focus on the effect of the slip velocity. Since the initial position of the 
cylinder is below the centerline, it is moving to the left in the lower region of the 
computational domain due to the given boundary condition (see Fig. [I]). For getting 
the slip velocity, we first compute the fluid horizontal speed on the streamline through 
the mass center in front of the cylinder at the distance of the half of the computational 
domain width and then minus the horizontal speed of the disk to obtain the slip 
velocity. The negative slip velocity means that the cylinder speed to the left is slower 
than the fluid speed to the left since the both signs are negative. We then have 
that the cylinder lags the fluid. For Re p =l, the slip velocity of the non rotating 



cylinder becomes negative after a short initial transition period as in Fig. 12 The 



one associated with free moving circular cylinder remains positive for, at least, the 
first 110 time units and then oscillates about zero. The cylinder with no rotation lags 
the fluid. But both cylinders migrate back to the middle between two walls due to 
that the wall repulsion force dominates the very weak slip velocity effect. We can see 
freely moving circular cylinder moves toward the centerline faster than the one with 



no rotation does in Fig. 12 When Re p =5, the slip velocity becomes negative after 



the initial transition period as in Fig. 12 The cylinder for the both cases lags the 



fluid. The slip velocity of the one with no rotation is about two times larger than the 
other one. The rotating velocity of each case of free motion is about constant speed 
after a short initial transition period. Fig. 13 shows that the relative velocity field to 
the horizontal velocity of the cylinder mass center in which we can clearly see that 
the cylinder lags the fluid. Both cases have stronger slip velocity which creates force 
pushing the cylinder to the region with faster flow speed, which is the region next 
to the bottom wall. Without the extra help from the Magnus lift, the one without 
rotation is closer to the wall. Hence when the initial position is not at the middle 
of two walls, the balance between the effect of the slip velocity, the Magnus lift and 
the wall repulsion does play a role for determining the equilibrium position of the 
cylinder. 



5. Conclusion 



We have investigated the motion of a neutrally buoyant cylinder of circular or el- 
liptic shape in two dimensional shear flow of a Newtonian fluid by direct numerical 
simulation. The numerical results are validated by comparisons with existing theo- 
retical, experimental, and numerical results, including a power law of the normalized 
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angular speed versus the particle Reynolds number. The rapid slow down of the 
normalized minimal angular speed of an elliptic cylinder is totally different from the 
behavior of the circular cylinder since the motion of an elliptic cylinder can transit 
from rotating to a fixed orientation when the particle Reynolds number is increased 
beyond the critical value. The midway between two walls is an expected equilibrium 
position of the cylinder mass center in shear flow. But when placing the particle away 
from the centerline initially, it migrates toward another equilibrium position between 
the wall and the centerline for higher Reynolds numbers which is caused by the in- 
terplay between the slip velocity, the Magnus force, and the wall repulsion force. The 
further study of the motion of a cylinder of axisymmetric shape in simple shear flow 
will be reported in a forthcoming paper. 
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